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Abstract. 

The Laser Interferometer Space Antenna will be able to detect the inspiral and 
merger of Super Massive Black Hole Binaries (SMBHBs) anywhere in the Universe. 
Standard matched filtering techniques can be used to detect and characterize these 
systems. Markov Chain Monte Carlo (MCMC) methods are ideally suited to this 
and other LISA data analysis problems as they are able to efficiently handle models 
with large dimensions. Here we compare the posterior parameter distributions derived 
by an MCMC algorithm with the distributions predicted by the Fisher information 
matrix. We find excellent agreement for the extrinsic parameters, while the Fisher 
matrix slightly overestimates errors in the intrinsic parameters. 
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1. Introduction 

Super Massive Black Hole Binaries (SMBHBs) are likely to be the most powerful 
sources of gravitational waves in the Universe. The Laser Interferometer Space Antenna 
(LISA) [Tj will be able to detect these systems out to the edge of the visible Universe. By 
studying SMBHB inspirals we may gain insight into the galaxy merger history and the 
role that black holes play in structure formation. SMBHB inspirals and their subsequent 
mergers also provide fertile ground for performing tests of general relativity j2] . 

To date, most data analysis development work for LISA has focused on Extreme 
Mass Ratio Inspirals j3] and the galactic confusion problem [H El El IZ|- SMBHB 
inspirals have received relatively little attention, perhaps in part because they are not 
expected to pose much of a challenge. The Markov Chain Monte Carlo (MCMC) method 
has emerged as a leading algorithm for LISA data analysis. The MCMC method 
can efficiently explore large parameter spaces, perform model comparisons, estimate 
instrument noise, while simultaneously providing error estimates for the recovered 
parameters. MCMC techniques have been applied to ground based gravitational wave 
data analysis '8\; a toy LISA problem jH]; and the extraction of multiple overlapping 
galactic binaries from simulated LISA data jS]. Here we make an initial foray into the 
SMBHB inspiral problem, and compare the parameter recovery accuracy of an MCMC 
search to the predictions of the Fisher information matrix. Considerable work has been 
done on using the Fisher information matrix to make predictions about LISA's resolving 
abihties for SMBHBs [121 HH lEl UHl IHj, so it is interesting to see how reliable those 
estimates might be. 

We find that the Fisher matrix approach yields very good estimates for the angular 
resolution and distance uncertainties, and that it tends to overestimate the errors in the 
component masses and the time of coalescence. 

2. The Model 

In this study we employ restricted post-Newtonian (PN) waveforms, which neglect 
the higher order harmonics, and treat the phase to 2-PN order ITH] . The 
detector response is modeled using the low frequency approximation ^01 ITB] . 
The gravitational waveform for a Supermassive black hole system consisting of 
two Schwarzschild black holes is described by a 9-D parameter set, x = 
{In(Mc), ln(/^), 9, 0, In(tc), l, fc, ^^{Dl),iP}, where Mc is the chirp mass, fj, is the reduced- 
mass, {6, (j)) give the sky location, is the time-to-coalescence, i is the inclination of 
the orbital plane of the binary, yj^ is the phase of the gravitational wave at coalescence, 
Dl'is the luminosity distance and if) is the polarization angle of the gravitational wave. 
We will describe the parameters D^, l, ipc, ^ as being extrinsic, while all the rest will be 
described as being intrinsic [Ej. We should mention that (6^,0) and tc would normally 
be classed as extrinsic, but they become quasi-intrinsic due to the motion of the LISA 
observatory. 
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Figure 1. A plot of the strain spectral density of the LISA response to a 10^— 10® Mq 
binary at z = 1. The dashed line indicates the RMS instrument noise level. 



We focus our attention on a particular SMBHB system consisting of a 10 — 10^ Mq 
binary system at z = 1. This gives corresponding values of {Mc, ^, Dl) = (4.93 x 
10^ Mq, 1.82 X 10^ Mq, 6.63 Gpc). The other parameters are defined by {9, 0, l, (pc, ^j) = 
(1.325,2.04,1.02,0.954,0.658) radians. We choose tc = 0.49 years, and set the time of 
observation to be Tots = 0.5 years. We used a sample cadence of 800 seconds. Due to 
the the fact that the equations describing the phase evolution of the wave break down 
before we reach the last stable circular orbit (LSO), we terminate the waveforms at 
R = 7M. The source has a signal-to-noise ratio of ~450. 

The one-sided noise spectral density for a LISA Michelson channel is given in the 
low frequency limit by 
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where L = 5 x 10^ kms is the arm-length for LISA, S'p^°^ = 4 x IQ'"^"^ / H z and 
gaccei — g X \{)~^^ rn? / s'^ / H z are the position and acceleration noise respectively. Using 
this formula, we generate instrumental noise from a Gaussian distribution. In Figure 
we display the power spectrum of the LISA response to the system we are investigating 
with instrumental noise from one channel of LISA. Our analysis uses LISA as a two 
channel detector, where the detectors are rotated by a factor of 7r/2 with respect to one 
another [101 . We also assume that LISA has a lower frequency cutoff of 10~^ Hz. 
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3. The MCMC approach 

In the Bayesian framework of data analysis, the goal is to calculate the posterior 
probability density functions for the parameters of the model, given the parameter priors 
and the data. However, this is extremely difficult, if not impossible, to do analytically. 
The advantage of the MCMC method is that it allows one to map out the posterior 
numerically. Our MCMC approach is implemented by using a Metropohs-Hastings 
algorithm to map out the posterior distributions of the 9-D parameter space for Super 
Massive black holes. The method works as follows: We assume that we have already run 
a search chain that has found the system we are looking for. We use the true parameter 
values as the starting point for our chain x. The algorithm then suggests a jump to a 
new point y using a proposal distribution q{'\x). We evaluate the Hastings ratio 
^ ^ 7^{y)p{s\y)q{x\y) 
7i(x)p(s\x)q{'y\x) ' 

and accept the jump with probability a = min(l,if), otherwise the chain stays at its 
current position. Here 7r{x) are the priors of the parameters and p{s\x) is the likelihood. 
We use uniform priors for all the parameters. The parameters In Mc, In In tc and 
In Dl were taken to be uniform in the range (—00, 00), cos 9 and cos l were taken to be 
uniform in the range [—1,1], and (pc were taken to be uniform in the range [0, 27r], and 
■0 was taken to be uniform in the range [0, tt]. The log- likelihood for a template h{x) 
given a signal s was assumed to have the form 

In p {s\x) — —-{s — h{x)\s — h{x)) . (3) 

Here the angular brackets define the standard noise weighted inner product. In order to 
achieve a healthy acceptance rate for the proposed jumps, we use a proposal distribution 
given by a multivariate normal distribution that is the product of independent normal 
distributions in each of the 9 eigendirection of the Fisher matrix, Vij{x). The Fisher 
information matrix describes the expectation value of the curvature of the log likelihood 
function at maximum likelihood: 



rjj(fML) = didj Inp ml) = ml) |9j^(^ml)) ■ (4) 

Here the over line indicates the expectation value. We used a generahzed notion of the 
Fisher matrix by employing the definition ^ijix) = {dih{x)\djh{x)) for points x away 
from maximum likelihood. The standard deviation in each eigendirection of I'ijix) was 
set to equal = l/sjDEi{x), where D = 9 and Ei{x) is the corresponding eigenvalue. 
The factor of 1/\/D ensures that the typical jumps are ~ la for the full multivariate 
normal distribution. In principle this choice of proposal distribution should yield a 
~ 69% acceptance rate for the proposed jumps. In practice we found an acceptance 
rate of ~ 33% for the system being studied. The lower acceptance rate is due to the 
Fisher matrix slightly overestimating the uncertainties in one of the eigendirections, and 
to a slight error in the estimate of the orientation of the corresponding eigendirection. 
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Figure 2. A plot of the MCMC chains for the 9 parameters describing the SMBHB 
binary during the initial burn-in phase. The straight line in each panel denotes the 
true parameter value. 



4. Results 

We started the chain with parameter values close to their true values, but offset by 
Ami/mi = -0.2%, Ams/ma = +0.3%, Az = +0.2, At^/tc = -2 x 10~^ A^ = -0.2, 
A0 = +0.2, At = +0.6, Alp = +0.2, and A0c = —0.4. The offset in each parameter 
was chosen to be ~tens of standard deviations away from the true source parameters 
(as estimated by the Fisher Information Matrix). In Figure |21 we see that the chain 
underwent a burn-in phase during the first 10,000 steps, before settling in close to the 
true source parameters. In repeated trials we found that the chain always settled in to 
the same region of parameter space after the burn-in was completed. When the chain 
was started off far from the true parameters the burn-in time became prohibitively 
long, so we do not recommend that the current algorithm be used for blind searches. 
We have developed a non-Markovian variant of the algorithm that can efficiently handle 
the search phase [IBj- The posterior distribution functions derived after the full search 
are identical to those found in the current paper. 

Following the burn-in, the chain was run for another 8 x 10^ points to explore the 
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Figure 3. A plot of the MCMC chains for the 9 parameters describing the SMBHB 
binary. The straight hne in each panel denotes the mean of the chain as calculated 
from the 8 x 10^ points. The dashed lines are the 2 — a predictions from the Fisher 
matrix evaluated at the mean of each chain. 



posterior distribution function. This took 8 days on a single 2 GHz processor. The 
first milhon points of the exploration phase of the chain are shown in Figure El We see 
that the MCMC chain moves around well, and show clear evidence of the high degree 
of correlation between the parameters Mc,fi,tc and ipc- The straight, solid line in each 
panel of Figure El denotes the mean value of the parameter, as calculated from the chain. 
Taking this mean value, we then calculated the standard deviation as predicted by the 
Fisher matrix at that point. This 1 — a error is given by 



(5) 

where Cij = (Fjj)"^. The two dashed lines in each panel denote the ±2a errors in each 
of the parameters. 

In Figure |31 we have plotted the 1-D marginalized histograms for each of the 
parameters. The solid line is the Fisher matrix prediction for the error at the mean 
of each chain. We can see that most of the extrinsic parameter histograms match 
the Fisher prediction almost perfectly. This is also the case for the {9, 0) parameters. 
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Figure 4. A plot of the marginalized 1-d histograms for each of the nine parameters. 
We see that most of the angular variables agree almost exactly with the prediction 
of the Fisher matrix. However, the remaining four parameters, which are all highly 
correlated, differ from the Fisher matrix prediction for the posterior. 



However, we can see that the histograms for (Mc, fi, tc, (pc) all differ from the prediction 
of the Fisher matrix. We attribute the error in the Fisher matrix prediction to the 
high degree of correlation between these parameters. In order to see just how correlated 
these parameters are, in Figure El we plot the 2-D marginalized histograms for the 
combinations {Mc,(pc), (a*, V'c), (tc'^c), and (McfJ,). We found that the Fisher matrix 
has a small eigenvalue in the (Mc, n, tc, (fc) direction which dominates the contribution 
to the proposed jumps in these parameters. The Fisher matrix tends to under estimate 
the eigenvalue in this direction, and it is also slightly off in predicting the corresponding 
eigendirection. Similar behavior was seen in all the examples we have looked at, which 
suggests that the Fisher matrix estimates for the uncertainties in the component masses 
found in the literature may be systematically high by about a factor of two. 

In a closely related study that appears in these same proceedings, Wickham, 
Stroeer & Vecchio ^se a Reversible Jump MCMC routine to study the posterior 
parameter distributions of a SMBHB system. They used a simplified model of the 
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Figure 5. This figure showing the 2-D marginahzed histograms demonstrates the 
large correlations between the parameters {Mc, fi,tc, fc)- 

SMBHB waveform, which removes the reduced mass from the parameter set. In contrast 
to our findings, they found that the Fisher matrix significantly underestimated the 
uncertainties in many of the model parameters. We do not know why our results are so 
different, but we do not think it is due to the differences in our MCMC implementations. 
The Reverse Jump method should yield the same results as our standard MCMC 
implementation when applied to models with fixed numbers of parameters. 

5. Discussion 

We have found that it is possible to construct a simple MCMC sampler for studying the 
posterior distribution functions for SMBHB inspirals in the LISA data streams. The 
marginalized posterior distributions recovered from the Markov chains are in good to fair 
agreement with the Fisher matrix predictions. In another work ^Hl, we have developed 
an advanced MCMC search routine that is capable of performing a blind search of the 
LISA data for SMBHB inspirals. 
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